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ABSTRACT 

Uniquely among the dwarf spheroidal (dSph) satellite galaxies of the Milky Way, Fornax 
hosts globular clusters. It remains a puzzle as to why dynamical friction has not yet dragged 
any of Fornax's five globular clusters to the centre, and also why there is no evidence that 
any similar star cluster has been in the past (for Fornax or any other dSph). We set up a suite 
of 2800 A^-body simulations that sample the full range of globular-cluster orbits and mass 
models consistent with all existing observational constraints for Fornax. In agreement with 
previous work, we find that if Fornax has a large dark-matter core then its globular clusters 
remain close to their currently observed locations for long times. Furthermore, we find pre- 
viously unreported behaviour for clusters that start inside the core region. These are pushed 
out of the core and gain orbital energy, a process we call 'dynamical buoyancy'. Thus a cored 
mass distribution in Fornax will naturally lead to a shell-like globular cluster distribution near 
the core radius, independent of the initial conditions. By contrast, CDM-type cusped mass 
distributions lead to the rapid infall of at least one cluster within At = l-2Gyr, except when 
picking unlikely initial conditions for the cluster orbits (~ 2% probability), and almost all 
clusters within At = lOGyr. Alternatively, if Fornax has only a weakly cusped mass distribu- 
tion, dynamical friction is much reduced. While over Af = lOGyr this still leads to the infall 
of 1-4 clusters from their present orbits, the infall of any cluster within At — 1-2 Gyr is much 
less likely (with probability 0-70%, depending on At and the strength of the cusp). Such a 
solution to the timing problem requires (in addition to a shallow dark-matter cusp) that in the 
past the globular clusters were somewhat further from Fornax than today; they most likely did 
not form within Fornax, but were accreted. 

Key words: stellar dynamics - methods: A^-body simulations - galaxies: kinematics and 
dynamics - galaxies: structure - galaxies: haloes 



1 INTRODUCTION 

The Fornax dwarf spheroidal (dSph) galaxy is the most massiv e 
undisrupted dSph satellite of the Milky Way dWalker et aI1l2009l) . 
Like all dSphs it is dark matter dominated even in its central re- 
gions. It is unique among the undisrupted dSphs in having globular 
clusters; it has five, with three of them at a projected distance out- 
side of the half light radius (see table [T}. There is also evidence of 
two shell-like structures, which may be the remnants of a merger 
occurring more than 2 Gyr ago dColeman et a l. 2004, 2 005b . 

One apparent paradox about these clusters is that, because they 
appear to orbit in a massive background of dark matter, they should 
be affected by dynamical friction which will cause their orbits 
to decay. Fornax's globular clusters are metal poor and very old, 
comparable to the oldest globular clusters in the Milky Way with 
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ages of the order of a Hubble time iBuonanno et al.l f 1998. 1999; 
Mac kev & Gilmord 12003d : iGreco et al .1120071) . During their life- 
time it would be expected t hat they fall to the centre of Fornax and 
form a nuclear star cluster dTremaine. Ostrik er & Spit zer Jr]|l975t 
lTremainelll976h . However, no bright stellar nucleus is observed in 
Fornax, or in fact any other dSph. This is known as the timing prob- 
lem for Fornax's clusters because it seems highly improbable that 
Fornax's globular clusters would be observed just briefly before 
they fall into the core. 

Several solut i ons to the timing problem have been proposed. 
lOh. Lin & Richer] j2000h suggested two ideas. First, that a pop- 
ulation of black holes transferred energy to the clusters through 
close encounters; and second, that a strong tidal interaction be- 
tween the Milky Way and Fornax could inject energy into their 
orbits. There is no observational evidence for a population of 
black holes in the centre of Fornax, while the currently ob- 
served proper motion indicates that the orbit of Fornax around 
the Milky Way never takes it closer than at present jDinescu et al.l 
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|2Q04| :[Lux. Read & Lake 2010). Thus, both ideas appear to be dis- 
favoured. lAngus & Diaferiol ( 120091) proposed that all but the most 
massive cluster could avoid sinking into the centre of Fornax if 
their current distance is much larger than projected but still within 
the tidal radius of Fornax of ~ 1.9kpc. However, this (i) requires 
special arrangement of the current projected positions and (ii) stills 
leaves a timing problem for the most massive cluster and is there- 
fore not a complete solution (moreover, their analysis was based on 
Chandrasekhar's simple dynamical-friction formula, which is not 
suitable for accurate estimates). 

Using numeri cal simulations and analytic arguments 
iGoerdt et alj J2006h proposed that the current distribution of 
the Fornax clusters can be explained by the diminution of 
dynamical friction on the edge of a cored matter distribution 
which causes the clusters to remain outside the dark matter 
core radius. Dynamical reasons fo r this ' core- st al ling' effect 
have been explored in | Read et all d2006h . Ilnoud d2009h . and 
ICole. Dehnen & Wilkinson! 201 1 ). Support for this result was pro- 
vided by [Sanchez-S alcedo. Reves-Iturbide & Hernandez! (2006) 
who showed that a cored matter distribution in dwarf galaxies can 
significantly delay the infall times of the globular clusters (even if 
Chandrasekhar's simple dynamical-friction formula is used). 

Measuring and/or constraining the dark-matter distribution 
in Fornax is interesting as a test of our current cosmolog- 
ical model. Collisionless cosmological simulations (which ig- 
nore the effects of baryons) predict a universal density distribu- 
tion for dark matter halos, with a central density cusp p oc r~ y 
where y~\ dDubinski & Carlberdfl99ll; In avarro, Frenk & White! 
[1991) . If the dark matter distribution in Fornax is found 
to deviate strongly from this prediction, this could imply 
that baryons have an important dynamical role in shap- 
ing the central dark matter d i stribution in dwarf galaxies 
(e.g. [Navarro. Eke & Frenkl Il996l |E1-Zant, Shlosman & Hoffman! 
iRead & Gilmord 120051. boerdt et all l2010l ICole etafl 
Pontzen & Governatoj 1201 2[) . or that we must turn to 
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more exotic cosmological models (e.g. Tremaine & Gunnj 19791 
Kochanek & Whitdl200(il.lHogan & Dalcantodl200d.lStrigari et all 
2006, Villaescusa-Navarro & Dalai 2011, Macci o' et al.ll2012l) . 



The first e vidence that dSphs may have a cons t ant de nsity 
core came from[klevna, Wilkinson, Gilmore & Evans (2003]) who 
found indirect evidence for a core in the Ursa Minor dSph. The 
Milky Way dSphs have been observed intensively in recent years, 
primarily because these systems are the most dark-matter domi- 
nated known. They contain mostly intermediate or old stellar pop- 
ulations which are likely to be well mixed in the dark matter po- 
tential because star formation ceased many dynamical times ago. 
This in turn implies that they are ideal laboratories for studying 
the mass structure of their dark matter halos. The intense obser- 
vational effort means that there is a wealth of kinematical data 
available to form the basis for theoretical models of these sys- 
tems. One line of approach has been based on the Jeans equations 
where a parametric light profile for the stars is assumed and a ve- 
locity dispersion profile is d erived based on a underlying param- 
eterised dark matter profile (IPenarrubia. Mc Connachi e & Navarro! 
l2008llStrigari et al.l2008l[Walker et al.l2009b. T his appr oach has its 
drawbacks (see for example lAmorisco & EvansluOllbl) . Most im- 
portantly, the degeneracy between the mass profile and the velocity 
anisotropy (which is poorly constrained), means that both cusped 
and cored density distributions are consistent with even the latest 
data. However, modelling the dSphs as two chemically distinct pop- 
ulations with dif ferent scale lengths, it appears that this degenera cy 
can be broken dBattaglia et al]|2008l : lAmorisco & Evanslboi lallbt 
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GC5 5.25±0.20 
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1.6 
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137 ±13 b -f 

138 ±8 s 
130.6±3.0 6 

136.1±3.1* -1.2±4.6 

135.5*3.1* 7.1 + 3.9 

134.0±6 c 5.9±3.4 

140.6±3.2'' 8.7±3.6 



Table 1. Data for the Fornax system. Column 3 shows the most likely mass 
from column 2. r c is the King (1962)model core radius for the globular 
clusters, and the half-light radius for the dSph. R p is the projected dis- 
tance of the cluster from the centre of Fomax. d\ os is the distance to each 
cluster and Av\ os the line-of-sight velocity relative to Fornax itself. Refer- 
ences: "Mackey & Gilmore (2003b),*Mackey & Gilmore (2003a), c Greco 
et al. (2007), rf Mateo et al. (1991), c Walker et al. (2009), / Buonanno et al. 
(1999),«Mateo (1998),''Buonanno et al. (1998). 



IWalker & Penarrubiall2oT l|); the results favour a cored density dis- 
tribution in the t wo dSph galaxies best analysed to date: Fornax and 
Sculptor (but see lBreddels et al.ll2012n . 

In this paper, we we follow the work of IGoerdt et al.l d2006l) 
by examining what the current location of the globular clusters can 
tell us about Fornax's mass distribution. Our work improves on this 
previous analysis in several key respects: (i) we use several mass 
models for the underlying potential in Fornax that sample the full 
range consistent with the latest data; (ii) we use the latest data for 
Fornax's globular clusters as constraints on their phase space dis- 
tribution; and (iii) we run thousands of /V-body models to sample 
the uncertainties in the cluster distribution and Fornax mass model. 

This large search of the available parameter space allows us to 
address whether or not there are multiple solutions to Fornax's tim- 
ing problem. To focus the discussion, we phrase the timing problem 
as a contradiction with either of the following two hypotheses 

(i) The Fornax globular cluster system is in a near-steady state, 
consequently none of the globular clusters should fall into the core 
of Fornax within a Hubble time. 

(ii) Our present cosmic epoch of observing Fornax is not spe- 
cial, consequently the system does not evolve significantly on a 
time scale short compared to a Hubble time. In particular, within 
1-2 Gyr none of the clusters should fall into the core of Fornax with 
high probability. 

The first hypothesis can be justified by the assumption that the 
present state of Fornax's globular cluster system is representative 
also for its past. This assumption, however, is not necessary, as one 
can easily think of alternative scenarios. The second hypothesis, on 
the other hand, is harder to avoid and similar to the cosmological 
principle, though here expressed in terms of the epoch of observa- 
tion rather than its vantage point and orientation. We will refer to 
contradictions with these hypotheses as the long-term and immedi- 
ate timing problem, respectively. 

This paper is organised as follows. Section|2]reviews the prop- 
erties of the Fomax system relevant for our study, Section |3]details 
our modelling approach, Sections |4]&(5]present the simulations re- 
sults and assess the probability that none of the five clusters will 
sink into the core of Fornax within either 1-2 Gyr or a Hubble time. 
Finally, in Section|6]we discuss the implications of our results and 
draw our conclusions. 
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Figure 1. Distribution of the Fornax globular clusters in mass and pro- 
jected distance from the centre of the Fornax dSph. The dashed vertical line 
indicates the stellar half-light radius of the dSph. 



in mass M and projected radius R v from the dSph. There are two 
interesting observations to be made from this figure. First, the radial 
distribution of clusters is consistent with that of the stars within 
Fornax: there is about half of the total cluster light within the stellar 
half-light radius of 668 pc. There are two possible interpretations 
of this. Either it is a coincidence, or the formation histories of the 
clusters and the Fornax galaxy are closely related, in particular, the 
clusters formed within the same entities as the stars. 

Second, there is a weak correlation between M and R p . In 
particular, the lightest cluster is furthest away from Fornax and 
the heaviest is the second closest. The remaining three are about 
equally massive and cover a spread of R p . This correlation is in the 
sense expected from mass segregation such as driven by dynamical 
friction. 



2 THE FORNAX SYSTEM 

We summarise in Table [T] the most relevant data for our study and 
their origin. 



2.1 The dSph 

As already discussed above, Fornax is the most massive of the 
Milky Way's dSph (except possibly for disrupted objects such as 
Sgr dwarf) and unique amongst (undisrupted) dSph to host a glob- 
ular cluster system. The very fact that Fornax can hold onto a 
globular-cluster system implies that Galactic tides cannot be strong 
enough to pull these clusters off. This in turn requires that its Galac- 
tic orbit never carries Fornax too close to the Milky Way, where the 
tidal forces become exce edingly strong. 

iLux. Read & Lake! fcoifjh estimate the peri-galactic radius of 
Fornax, based on its observed position, distance, radial velocity and 
proper motion, to be 100-130 kpc, which indeed is only slightly 
smaller than its current distance of about 140 kpc (see TableQ}. 

Based on this result, we estimat e the tidal radius for Fornax 
using the method o f lReadetaljj2006h , where the dSph is modelled 
as a spherical sa tellite orbiting the Milky Way repres ented by a 
Hernquist ( 1990) model. We then solve equation (7) of iRead et al.l 
(2006), which accounts for the orbit of the satellite about the host, 
and the orbits of the stars within the satellite. We find a tidal radius 
of 1.8-2.8 kpc, based on a the range of masses for Fornax given in 
Tableland using the extremal values for the the orbital data taken 
from iLuxet al. (2010) and a total (extended) mass for the Milky 
Way of l-2xl0 12 M o . 



2.2 The globular clusters 

O ur principle sources for glob ular-cl uster data are those published 
bv lMackev & Gilmord d2003al lbh and lGreco et al.l d2007h who have 
carried out thorough surveys of the Fornax globular clusters. For 
our purposes the main data required are the clusters' masses, sizes, 
three dimensional positions and velocities. 

The best estimates for these quantities are given in Table [T] 
The values for the core radius r c of eac h cluster are based on the sur - 
f ace bright ness profiles calculate d in iMackev & Gilmorel d2003bh. 
These are lElsonTF all & Freeman (1987, EEF) models and the Kind 
d 19621) model core radius r c is related to the EFF scale parameter a 
by r c = a(2 2/y - 1) 1/2 , where y is the power law slope of the surface 
brightness at large radii. 

Fig.[T]plots the distribution of the five Fornax globular clusters 



3 MODELLING APPROACH 

The basis for our approach is to take the most up-to-date observa- 
tions of Fornax's globular clusters and combine these with plausi- 
ble mass models consistent with the latest kinematic data for For- 
nax's stars. We then create, for each Fornax mass model, initial 
conditions for the Fornax globular-cluster system, which are con- 
sistent with the relevant observations, and evolve them for ten Gyr 
into the future. 



3.1 Mass models for Fornax 



IWilkinson etafl d2002h and lKlevna etafl 120021) demonstrated that 
the mass-anisotropy degeneracy inherent in kinematic modelling 
can be broken using distribution-function modelling of sufficiently 
large kinematic data sets. To take full advantage of the recent data 
set of more than two t housand individual stellar velocities in For- 
nax dWalker et al . 2009), Wilkinson et al. (in preparation) apply the 
Markov-Chain-Monte-Carlo (MCMC) technique to dynamical and 
mass models of Fornax. The mass profile and stellar-luminosity 
profile are modelled independently using spherical double-power- 
law profiles of the form 



/Wk-i(''> = a> I ■H ( 1 + (t; 



(1) 



The stellar distribution functions are ca lcula ted numerica lly fol- 
lowing the approach of iGerhardl dl99ll) and ISahal d 19921) allow- 
ing various velocity anisotropy profiles. A s in the earlier work 
(Wilk inson et"ai]|2002l : iKlevna et alj [2002). the models are com- 
pared to the data on a star-by-star basis. 

Further details of the modelling and its results will be pre- 
sented elsewhere. Here, we simply use the above MCMC model 
ensemble to inform our choice of halo models. Rather than consid- 
ering a single best-fit model of Fornax, we select four mass models, 
each of the form |[T) but truncated at very large radii via 

p(r)=p mM (r) sech(r/10kpc), (2) 

and with parameter values as detailed in Table [2] These models 
span the range of models consistent with the kinematic data. The 
three models WC (weak cusp), IC (intermediate cusp) and SC 
(steep cusp) have parameters y Q , y m , 77, r s , and Af«, directly taken 
from the MCMC chain outputs, and refer to the highest likelihood 
models with density slope 

dlnp 



y(r) = - 



din j- 



(3) 
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Model 


Name 


70 


Too 


V 


M„ 


'"s 


yioopc 


M(1.8kpc) 


LC 


large core 


0.07 


4.65 


3.7 


8.00x10 s 


1.4 


0.1 


4.12x10 s 


WC 


weak cusp 


0.08 


4.65 


2.77 


1.23 x10 s 


0.62 


0.1 


1.03 x10 s 


IC 


intermediate cusp 


0.13 


4.24 


1.37 


1.51x10 s 


0.55 


0.5 


1.03 x10 s 


sc 


steep cusp 


0.52 


4.27 


0.93 


1.98x10 s 


0.80 


1.0 


1.07 x10 s 



Table 2. Parameters for halo mass models (equation [T] used in the simulations (yo, y m , r], Moo, r s ), 
as well as the resulting logarithmic density slope y(r) = -dlnp/dlnr at lOOpc and the mass within 
1 .8kpc, our lower limit for the tidal radius of Fornax. Radii are given in kpc and masses in M G . 




R [kpc] 

Figure 2. Radial profiles of density (top), enclosed mass (middle), and 
logarithmic density slope (bottom) for the four mass models used in our 
simulations (see also Table|2j.The data points the middle panel correspond 
to the mass estimates by Walker & Penarrubia (2011)for two chemically 
distinct sub-populations. 



of, respectively, 0.1, 0.5, and 1.0 at r= lOOpc. 

The fo urth model, LC (large core), was motivated by the re- 
cent work of Walker & Penarrubia (2011). These authors applied a 
non-parametric statistical modelling technique to two chemically 
distinct stellar populations within Fornax to define the enclosed 
mass at the half light radii of the two populations. The resulting 
model possesses a large core with near-constant density and y x 0. 1 
forr<500pc. 

The radial profiles of density, enclosed mass, and logarith- 
mic density slope of the four mass models are shown in Fig. [2] 
for comparison. These models cover a wide range of inner den- 
sity slopes, includin g shallow profiles, such as suggested by 
(Gilmore et al. 2007) based on observations of dSph galaxies, but 
also a steep cusp, such as predicted by cosmolo gical simulations 
dPubinski & Carlberdfl99ll iNavarro et al.lll996l) . Note, however, 
that our models represent the overall mass distribution of Fornax 
including both the stars and the dark matter. 

Figure [3] shows the stellar velocity dispersion for of our three 
MCMC-based Fornax models (curves) pl otted together with t he ob- 
served velocity dispersion as measured by Walker et al. (2007|). The 
model velocity dispersions were calc u lated assuming an ergodic 
distribution function with a iPlummerl d 1 9 1 If) density profile with 
core radius 668 pc for the stellar component, but the respective halo 




weak cusp (WC) 

intermediate cusp (IC) 

steep cusp (SC) 



0.5 1 1.5 2 

r [kpc] 

Figure 3. The observed stellar velocity dispersion for Fornax as a func- 
tion of projected radius (Walker et al. 2007), and the simple predictions ob- 
tained for the four mass models used in this study under the assumption 
of isotropic velocities. The purpose of this comparison is not to assess the 
relative merit of the mass models, but merely to demonstrate that their nor- 
malisations are reasonable. 

model for the underlying mass distribution. In view of the fact that 
these simple ergodic models have not been fitted to the velocity- 
dispersion data (apart from the assumed mass models), they pro- 
vide surprisingly good a description of these data. (This and the 
similarity between the model predictions is exactly the reason why 
inferring the mass profile from data like these is hardly possible.) 

The large core model (LC) has been normalised to roughly 
agree with the estimates for the enclosed mass derived from two 
different tracer populations by Walker & Penar rubia] fcoill) and 
plotted in Fig. [2] 

3.2 Modelling the globular cluster system 

If we compare the distance to the Fornax dSph with the distances 
to the individual clusters in Table [TJ it can be seen that the mea- 
surements of the distances are not accurate enough to provide re- 
liable three-dimensional locations within the dSph. We therefore 
draw for each simulation random line-of-sight distance offsets to 
Fornax from a uniform distribution between and 2 kpc, the ap - 
proximate tidal radius of the system JWalker & Penarrubiall201 II) . 
For the velocities, we choose a similar statistical approach by sam- 
pling the full space velocity from a bi-variate Gaussian distribution, 
specified by the total velocity dispersion cr and the anisotropy pa- 
rameter 



This is done in such a way that for clusters GC2-5 the line-of-sight 
velocity matches the observed value (which may be considered a 
prior for our sampling). 
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As the spatial distribution of the clusters is consistent with that 
of the stars in Fornax, it seems reasonable to base our kinematical 
parameters for the clusters on the observed stellar kinematics. The 
measured stellar velocit y dispersion is approximately flat over the 
range of radii observed dWalker et al.ll2007t lLokasll2009l) . and we 
use cr = 10.5 km s - '. For the v elocity aniso tropy we assume /3 = 
-0.33 as suggested for the stars dLokasl2009h , i.e. a mild tangential 
bias, which gives cr, as 9.5km s _I and cr, a; 1 1 km s _I . 

With this modelling approach, in particular the wide range of 
initial radii sampled, we generate many different orbits for the glob- 
ular clusters, covering the complete range of all possible orbits for 
them. By allowing such a wide distribution of cluster orbits, we 
can explore the effects of possible narrow choices for the cluster 
distribution function, such as preferably circular orbits , afterwards 
by restricting our analysis correspondingly. 



4 RAW SIMULATION RESULTS 



As detailed at the end of 93.21 our simulations cover a wide range 
of initial cluster orbits, some of which may not be very realistic. 
In this section we ignore any implications of our sampling of the 
initial cluster orbits and simply consider the individual simulations 
on their own merit. 

For each simulation, we need to quantify in how much each 
simulated cluster has suffered from dynamical friction and has 
sunken into the core of the dSph. In this respect it is clearly bet- 
ter to use an orbital property instead of an instantaneous quantity, 
such as radius (projected or intrinsic). Therefore, we use as our 
main characteristic the apo-centric radius r apo of the instantaneous 
cluster orbit. This is obtained from the cluster's instantaneous an- 
gular momentum L and specific energy E as the larger of the two 
radii for which 



E= — +<SKfy 
2r 



(7) 



3.3 A'-body simulations 

For each Fornax mass model, we ran 700 A-body simulations, each 
with different initial conditions for the five clusters. 

The individual cluster positions and velocities are drawn as 
described in the previous sub-section and their masses are taken 
from Table [TJ The clusters are represented by individual massive 
but softened particles with density profiles 



p(r)oc(r +e ) 



-7/2 



(5) 



where e = 5 pc, comparable to the cluster core radii. 

To generate the A'-body initial conditions for the Fornax mass 
models, we sample positions from the density ([2j and veloci- 
ties from self-consistent ergodic distribution functions /(£), which 
only depend on the specific orbital energy E, thus giving every- 
where isotropic velocity distributions. The forces between particles 
representing Fornax are softened with softening length e= lOpc. 

We enhance the resolution of the TV-body model in the inner 
parts (where dynamical friction occurs) by increasing the sampling 
probability by a factor g(E)~ l which is compensated by setting par- 
ticle masses fij proportional to g(Ej). We used 



8(E) oc 



l+qr". (E) 

1 cue v 

A (E) + r'! 

are / s 



(6) 



with q = 4 the ratio between maximum and minimum particle mass 
and ^circCE) the radius of the circular orbit with specific energy E. 
Testing this method for our particular purposes, we found that it al- 
lows a reduction of N to half at the same central resolution without 
any adverse effects. Based on convergence tests of decaying cluster 
orbits (see Appendix [A}, we use N = 2x 10 6 particles to sample 
each of the Fornax mass models. 

The mass ratio of the lightest background particle in all our 
models to the lightest cluster is 1:1680 (in model WC) and that 
of the heaviest particle to the lightest cluster 1:62 (in model LC). 
The clusters orbit mainly within the high-resolution region of our 
models defined by the volume where the resolution is better than 
would have been achieved with the same number of particles of 
constant mass (within approximately 2kpc). 

The simulations are performed using th e publicly av ailable A'- 
body code gyrfalcON, which uses Dehnen's J2000ll2002h 0(N) al- 
gorithm for force approximation. The total energy conservation was 
typically a few parts in 10 4 . 



(the smaller is the peri-centric radius). O(r) is the instantaneous 
gravitational potential of the A-body model, estimated using spher- 
ical averaging. For a small fraction of simulated clusters, the initial 
orbits were unbound. Such simulations are discarded only for anal- 
ysis of the affected cluster orbit, but not for the others. 

For each of the four mass models of Table [2] Fig. |4]plots the 
distributions of each cluster orbit in initial and final r apo at / = 2Gyr 
and t = lOGyr in the left and right sub-panels, respectively. 

4.1 Orbital decay after 2 Gyr 

After 2 Gyr (left sub-panels in Fig. [4} r apo is significantly less than 
initially for most orbits with initial r apo < 2- 3kpc, except for model 
LC, which we discuss separately in 94.31 

However, both GC 1 (red) and GC5 (cyan) rarely sink substan- 
tially, and GC 1 only ever falls into the centre of Fornax for the most 
cusped model SC. This is not surprising since GC1 is not only by 
far the lightest of the five clusters (see Table [T}, but also the one 
furthest away (in projection) from the centre of Fornax and hence 
less likely to pass through high density regions. As such it requires 
most dynamical friction to fall into Fornax, but will suffer least, 
since drag force oc mass 2 p. Though GC5 is five times more mas- 
sive, it has the second greatest projected distance from the centre 
of Fornax and hence also suffers significantly less friction than the 
other clusters for most orbits sampled. 

The remaining three clusters GC2-4 are all dragged inwards 
when initially r apo < 2 - 3 kpc, presenting the Fornax timing prob- 
lem. For all of these clusters, orbits with initial r apo < lkpc show 
similar effects of dynamical friction, presumably because orbits 
with apo-centres as small as that spend sufficient time in high- 
density regions to suffer substantially from dynamical friction. 

Since R p < r apo , an initial r apo ~ lkpc is the absolute possi- 
ble minimum for GC2, which is observed at R p = 1.05 kpc. Conse- 
quently GC2 only rarely falls in as much as GC3 and GC4. 

Apart from the initial r apo , the infall of these clusters depends 
most strongly on the inner density profile of the background mass 
distribution. Model SC shows the greatest effect of dynamical fric- 
tion on the clusters. For GC3 and GC4, a significant proportion of 
the simulations find their instantaneous apo-centres inside 30 pc. 
GC2 does not show such a marked effect but a number of simula- 
tions already have an apo-centre inside 100 pc. 

By contrast, model WC shows significantly reduced dynam- 
ical friction. GC3 again is most affected. However, even in the 
most extreme cases, the apo-centres have decayed to no less than 
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Figure 4. Apo-centric radii r ap0 of the instantaneous cluster orbits after t = 2 and 10 Gyr for all simulations per halo model (as indicated) plotted versus r apo 
of the initial orbit. The shaded region indicates the most likely value for the current tidal radius of Fornax. Any initial r apo greater than that would have been 
strongly affected (likely removed from Fornax) by the Galactic tidal field (not modelled in our simulations). The thin horizontal lines indicate the observed 
projected radius S D b s for each cluster, which is necessarily less than its initial r apo . Any final r apo below these lines would be inconsistent with the present 
cluster location. 



~ lOOpc from the centre. In this model the logarithmic density 
slope y (equation^ was initially only 0.1 at 100 pc, which implies 
a near-flat density profile within this radius. As expected, model IC 
is intermediate between models WC and SC. 

4.2 Orbital decay after 10 Gyr 

After 10 Gyr, the trends shown at 2 Gyr are amplified. Four of the 
globular clusters fall into Fornax for most simulations. Only GC1 
still shows not much evidence of migrating to the centre of Fornax. 
All cluster have a small fraction of simulations where they remain 
at large r apo . This requires the initial r apo to be large as well. In gen- 
eral, r apo decreases at least slightly, even when initially large (the 
few significant increases of r apo are caused by cluster encounters). 

For the most cusped model SC and to some degree for model 
IC too, GC3, GC4 and even GC2 can obtain values for r apo down 



to close to their softening length of 5 pc — i.e. they are essentially 
at the very centre of the galaxy. For model WC, this is not the case, 
i.e. none of the clusters can reach the very centre of the galaxy in 
this case. However, the vast majority of simulations do not obtain 
such very small r apo . 

We observe that the distributions of r apo after 10 Gyr are quite 
similar between models WC, IC, and SC (in particular the latter 
two), when GC3 has settled at r apo < lOOpc for most of our sim- 
ulations, in fact all simulations with initial r apo < 2.8 kpc, and the 
distributions for r apo of GC4 are remarkably similar. For the most 
massive GC3, this similarity reflects the fact that the cluster has 
essentially sunken to the very centre of the galaxy. 

However, this cannot explain the similarities for GC4, which 
we think is caused by a similarity of the background mass pro- 
files after 10 Gyr. The action of cluster dynamical friction causes 
dynamical heating of the background particles dragging the clus- 
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Figure 5. Ths initial (black) and final (t = lOGyr, coloured) density profile 
of six halo models (stacked /V-body models) from simulations where either 
three clusters sunk into Fornax or none. The dashed curves represent ±lcr. 



ters. This in turn erase s the initial central de nsity cusp dRead et all 
l2006tlCole et alj|201 lh . lGoerdt et at] d2010h discuss the formation 
of a density core in this way and provide an empirical formula for 
the expected core size created as the clusters fall to the centre. For 
model WC this predicts a core radius of 262 pc which is larger than 
our stalling radii but agrees within a factor of a few. However this 
formula fails to predict the behaviour of the clusters in the IC and 
SC models as it gives stalling radii of 176 pc and 100 pc respec- 
tively. 

In Fig. [3] we plot the density profiles of the final models for 
cases with and without cluster infall. For the steep and intermediate 
cusp models, the central density profiles are significantly reduced 
and in fact have become much more similar to each other. However, 
only for model SC is this reduction slightly stronger when clusters 
have reached the core of Fornax. In all cases, the density still keeps 
on increasing inwards, and hence does not really correspond to a 
constant-density core in the classical sensfl 

4.3 The large-core model 

The large core model (LC) shows very unusual behaviour. After 
2Gyr, all globular clusters have apo-centres closely clustered to- 
gether with a strong peak at ~ 1 kpc. This becomes even more pro- 
nounced at lOGyr. Detailed examination of the cluster orbits in 
individual simulations shows two interesting behaviours. First, or- 
bits with initial r apo > 900pc decay move in quite rapidly (mostly in 
less tha n 2 Gyr) to r ap o ~ 900 pc, where they stall. T his confirms the 
work of iGoerdt etafl (2006) and [Read etaTI j2006l) which showed 
that massive satellites orbiting outside of an harmonic density core 
stall at the edge of the core. This behaviour is believed to be due 
to the reduction of dynamical friction due to the resonant effects of 
particles in the harmonic core. 

1 For model WC, the density actually increases slightly at r < lOOpc. 
This is presumably caused by a slight instability of this model, which has 
d// dE > for some range of sp ecific energies E and hence may be unstable 
(see lBinnev & TremainehoOSi §5.5). 
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Figure 6. The evolution of the orbit of GC4 in one of our simulations with 
the large core (LC) model. The top and middle panels show the development 
of the orbital radius and specific energy for GC4, respectively. The latter 
grows after approximately 2Gyr as the orbit moves outwards. The overall 
conservation of energy (bottom panel) is not affected by this orbital change. 



Second, any cluster which has an initial orbit within the har- 
monic core move out to the edge of the core. Fig.|6]shows the evo- 
lution of the orbit of GC4 in one of our simulations. As can be seen 
the cluster orbit absorbs energy as it expands. The overall conserva- 
tion of energy for the simulation is not affected by this change and 
energy is conserved to approximately 3 parts in 10 4 during the sim- 
ulation. This behaviour is unexpected and we do not believe that it 
has been reported previously, though there is some evidence for or- 
bital radii expanding again after falling in at the edge of harmonic 
cores in what has been called the 'kickback effect' dGoerdt et all 
l201Ct l Inoue 2009). We will discuss this further in section|6] 



4.4 Summary 

The raw empirical results from our simulations can be summarised 
as follows. 

(i) Cluster orbits with large initial r apo are not significantly af- 
fected by dynamical friction, because these orbits spend no or too 
little time in high-density regions, where frictional drag is exerted. 
However, most simulations which initial r apo less then the current 
tidal radius of Fornax suffer from dynamical friction. 

(ii) Cluster GC3 is most likely affected by dynamical friction, 
followed by GC4 and GC2, while clusters GC1 and GC5 are 
least likely affected after 2 Gyr. This ordering is expected from the 
masses and initial projected radii of the clusters. 

(iii) For all except model LC (large core), cluster GC3 always 
reaches the core of Fornax within 10 Gyr (unless its initial orbit 
was unrealistically large with r apo > 2. 8 kpc beyond the present-day 
tidal radius of Fornax), constituting the long-term timing problem. 

(iv) The effect of dynamical friction at 2 Gyr is increasing with 
the central mass density from model WC to SC, as expected from 
Chandrasekhar's dynamical friction formula. 

(v) The effect of dynamical friction after 10 Gyr is more similar 
for the three halo models with weak to steep cusps than after 2 Gyr. 
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This similarity can be understood, at least qualitatively, by stalling 
of dynamical frictions as consequence of core formation. 

(vi) Model LC shows no effect of dynamical friction, but 
rather the opposite: 'dynamical buoyancy', when the clusters are 
pushed out of the core. Like dynamical friction, this effect appears 
strongest for the most massive cluster. 

In particular, for a CDM-type steep cusp (model SC), almost all 
simulations with initial r ap0 < 2.8kpc (Fornax's tidal radius) suffer 
significantly from dynamical friction within 10 Gyr, or even within 
only 2 Gyr. This reflects t he Fornax timing problem and is in contra- 
diction with the claims o fl Angus & Diaferio|j2009l) (as discussed in 
the introduction), possibly because our sampling discourages very 
nearly circular orbits, but the usage of Chandrasekhar's simple for- 
mula by these authors certainly plays a role too. 



5 THE PROBABILITY OF CLUSTER SINKING 

Our results presented in the previous Section and Fig. |4] demon- 
strate that the more massive clusters GC3 and GC4 will fall into 
the core of Fornax within 2 Gyr or less, unless either the mass dis- 
tribution of Fornax has a core or the cluster has an initial orbit with 
large apo-centric radius r apo . This latter scenario requires that the 
observed cluster position is either near peri-centre (special orbital 
phase) or has intrinsic radius much larger than projected (special 
projection geometry). Either is atypical, implying that this scenario 
is inherently unlikely. In this section we try to quantify just how 
unlikely. 

Since we know neither the current orbital phase nor projection 
geometry for any of the globular clusters, our most sensible ap- 
proach is to marginalise over both, assuming uniform distributions. 
However, as we will see, our sampling of initial cluster position 
and velocity did not generate a uniform sampling in these quanti- 
ties. As a consequence, any quantitative statements, based on the 
raw results, about the probability of cluster infall will be biased. In 
order to eliminate this bias, we must emulate a uniform sampling 
of initial orbital phases and projection geometries. 

To this end, we need a quantity for each simulated cluster 
which would follow a known distribution were orbital phase and 
projection angle drawn randomly. A natural such quantity is the 
fraction 

p(fl<i?p|orbit) 

of orbital phases and projection geometries for which the projected 
radius R is smaller than actually observed for the corresponding ini- 
tial cluster orbit (see Appendix |B]f or a formula and its derivation). 
Under our basic assumption of random orbital phase and projec- 
tion, p(R < i? p |orbit) is uniformly distributed between and 1, in 
particular p(R < 0|orbit) = and p(R < r apo |orbit) = 1. We can there- 
fore emulate a uniform distribution in orbital phase and projection 
geometry by a uniform distribution in p(R < R p |orbitf|. 

In Fig. [7j we plot for each simulated cluster and each halo 
mass model the instantaneous apo-centric radius r apo after 2 and 
10 Gyr against p(R < _R p |orbit). For all mass models, there is a clear 

2 In addition to p(R < R p |orbit) one may also utilise the fraction p(v z < 
Ml os |i? p , orbit) of the line-of-sight velocity to be less than observed. This frac- 
tion should also be uniformly distributed independently of p(R < R p |orbit), 
and thus allow an additional independent constraint. However, its computa- 
tion is more involved and it seems unlikely that much would be gained from 
it, mainly because of the relatively large uncertainty of the observed v\ os . 



correlation between p(R < /? p |orbit) and r apo at later times in the 
sense that larger r apo (r > 0) are achieved only when the initial pro- 
jected radius was relatively small for the initial orbit. This makes 
perfect sense: for small p(R < /f p |orbit) the orbit spends most of its 
time at r > R p with little dynamical friction. 

What is somewhat surprising, however, is how strong the cor- 
relation between p(R < W p |orbit) and r apo (r > 0) actually is, given 
that the initial orbits cover a wide range of eccentricities. For sim- 
ulated clusters with initial r apo < 2.8kpc, we have split their initial 
orbits into low and high eccentricity 



'"apo '"peri 

with open symbols and crosses in Fig.UJcorresponding to e < 0.4 
and e > 0.4, respectively. For the halo models IC and SC, we can 
see some differentiation between these two groups of initial orbits, 
in particular at t = 2 Gyr, in the sense one would expect: eccentric 
orbits obtain smaller r apo (because they have smaller initial r peri and 
hence suffer more dynamical friction) unless the observed R was 
initially untypically small (when they spend most of their time at 
large radii). 

As is evident from Fig. [7] the distributions of p(R <_R p |orbit) 
from our simulations are not uniform: there are many more sim- 
ulations with small p(R < i? p |orbit), in particular for clusters with 
small R p , such as GC3 and GC^E This non-uniformity is simply a 
consequence of our sampling procedure, which favours orbits with 
r nfo ^ Rp (in particular for clusters with small ^ p , such as GC3 and 
GC4) resulting in non-uniform sampling of orbital phase and pro- 
jection geometry and hence introducing a bias. 

Since r apo (f > 0) is strongly correlated with p(R < i? p |orbit), we 
can read off Fig. |7j the unbiased probabilities for infall into For- 
nax. For example, in model SC r apo (? = 2Gyr) > 200pc for GC4 
(magenta) in almost all simulations with p(R < i? p |orbit) < 0. 1 but 
hardly for any simulation at larger p(R < 7? p |orbit), implying that 
cluster GC4 will have r apo > 200 pc after 2Gyr with ~ 10% proba- 
bility in halo model SC. 

A more quantitative procedure for removing the bias of our 
non-uniform sampling of orbital phase and projection is to give 
each simulated cluster orbit a weight such that the weighed distri- 
butions of simulated orbits have, or at least are consistent with, a 
uniform sampling. Obviously, the weight to choose is simply the 
reciproce of the actual frequency f(p(R < 7? p |orbit)) of fractions 
p(R < K p |orbit) for all simulations of the same cluster in a particu- 
lar halo model. The probability for, say r apo (f = 2Gyr) > 200pc for 
GC4, is then obtained as the weighted fraction of simulations that 
obtain r apo (r = 2 Gyr) > 200 pc for this cluster. 

Combining this for all clusters, we estimate the probability 
that none of the five cluster when initially on orbit with r apo < 
2.8kpc (our upper limit for the tidal radius of Fornax) will fall into 
Fornax, in the sense that r apo < lOOpc or < 200pc, within 1, 2, or 
10 Gyr. The results for each halo model are given in Table ^(ex- 
cept for model LC when no cluster ever falls into Fornax). Columns 
2-4 give the probabilities based on all of our simulations. If we re- 
strict the analysis to high or low eccentricities (columns 5-7 and 
8-10, respectively) or to orbits with initial r apo < 1.8kpc (not given 
in Tabled, the results are hardly altered. 

In other words: these probabilities depend mainly on the halo 

3 For model LC, no simulation has small p(R < S p |orbit) for clusters with 
R p > 1 kpc. This is a consequence of the larger mass of this model which 
implies that our sampling did not obtain any nearly unbound orbits with 
large initial r ap0 . 
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Figure 7. Instantaneous apo-centric radius r ap0 at t = 2 and lOGyr (as in Fig.[4j for all simulated clusters in a given halo model (as indicated) plotted versus 
the fraction p(R < S p |orbit) of orbital phases and projections for which R^R p on the initial cluster orbit. A uniform distribution of p(R < R p |orbit) e [0, 1] 
corresponds to an unbiased distribution of orbital phases and projection angles. Dots are for orbits with initial r apo > 2.8kpc, crosses for eccentricity e > 0.4 
and circles for e < 0.4. The horizontal lines indicate R p and the shaded band indicates the estimates for Fornax's tidal radius. 



model and only weakly on the distribution function of the cluster 
orbits. This is a direct consequence of the tight relation, seen in 
Fig. between p(R < /f p |orbit) and r apo at later times. This result 
implies that we do not need to investigate further the implications 
of different distributions functions for the initial cluster orbits. 



6 CONCLUSIONS 

The Fornax galaxy is unique among the Milky Way dSphs in having 
five surviving globular clusters. These clusters are metal poor and 
very old - comparable with the oldest globular clusters in the Milky 
Way. There is a well-known timing problem for these clusters: over 
a Hubble time or even a much shorter time scale dynamical friction 
should drag them to the centre of Fornax, where they would form 
a nuclear star cluster. Yet no such stellar nucleus is observed in 
Fornax or, in fact, in any other dSph galaxy. 



In this study, we have extended previous work on what the 
current location of Fornax's globular clusters can tell us about For- 
nax's mass distribution. We explored four different mass models 
for the underlying potential in Fornax; we used the latest data for 
Fornax's globular clusters as constraints on their phase space distri- 
bution; and we ran thousands of ,/V-body models to sample the un- 
certainties in the cluster position and velocity distribution for each 
of our five mass models. This large search of the available parame- 
ter space allowed us to hunt for viable solutions to Fornax's timing 
problem. 

6.1 Caveats 

Our models all assume a spherical mass distribution for Fornax, 
yet the stars are observed to be ellipsoidal in projection, implying 
that the true intrinsic distribution is aspherical, most likely triaxial. 
The parameter space of such configurations is considerably larger 
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any e e < 0.4 e > 0.4 

model 1 Gyr 2 Gyr 10 Gyr 1 Gyr 2 Gyr 10 Gyr 1 Gyr 2 Gyr 10 Gyr 







probability for no cluster with r apo 


< lOOpc 


WC 
IC 

sc 


1 

0.17 
0.10 


1 <0.001 1 1 <0.001 
0.032 <0.001 0.15 0.022 <0.001 
0.017 <0.001 0.082 0.014 <0.001 


1 1 <0.001 
0.19 0.037 <0.001 
0.11 0.020 <0.001 






probability for no cluster with r apo 


<200pc 


wc 

IC 

sc 


0.92 0.32 <0.001 0.87 0.29 <0.001 
0.087 0.015 <0.001 0.049 0.01 1 <0.001 
0.076 0.014 <0.001 0.063 0.01 1 <0.001 


0.97 0.37 <0.001 
0.11 0.018 <0.001 
0.086 0.015 <0.001 



Table 3. Probability, estimated as the fraction of orbital phases and projec- 
tion geometries (consistent with the observed projected cluster positions) 
that no cluster with initial r apo < 2.8kpc (Fornax's tidal radius) has reached 
'apo < lOOpc or 200pc after 1, 2, or 10 Gyr of simulation, depending on the 
mass model and the eccentricity of the initial orbit. 



(and the freedom in cluster orbital projections smaller), potentially 
allowing more solutions to the timing problem. However, previous 
studies suggest that, while dynamical fr iction is on average stronger 
on b ox orbits than on loop orbits dCapuzzo-Dolcetta & Vicaril 
2005), triaxiality h as no sifnificant overall effect on the strength o f 
dynamical friction l lCora. Vergne & Muzzio||200ll:lsachaniall2009h . 
This can be understood as cancellation of two opposing effects. 
First, because angular momentum is not conserved, there is no bar- 
rier for eccentric orbits to reach arbitrary small radii and hence high 
densities and strong drag. Second, also because angular momentum 
is not conserved, successive peri-centric passages have different ra- 
dius and density, such that an orbiting cluster only rarely suffers 
very high drags. 

Another minor caveat of our modelling is our ignorance of the 
tidal field of the Milky-Way potential. In principle, it would be triv- 
ial to have our A'-body models orbit in a static model for the Galac- 
tic potential. However, such a procedure would add even more only 
weakly constrained parameters without adding significant benefits. 
With our existing models, we take the Galactic tidal effects into ac- 
count by discounting any simulated cluster orbits with apo-centric 
radius beyond (best estimates of) Fornax's tidal radius. 

Finally, by modelling each globular cluster as a massive parti- 
cle, we have ignored their inner dynamics and tidal interaction with 
Fornax. Mass-loss rates for these clusters, how ever, are likely to b e 
too small to significantly reduce their mass dGoerdt etal.1 (20101) . 
Tidal disruption is arguably only relevant for GC1, as discussed in 
aObelow. 



6.2 Solutions to the timing problem 

Our simulations demonstrate that for normal mass models (WC, 
IC, or SC) for Fornax the infall of clusters GC3 and GC4 within 
a Hubble time is unavoidable and the infall of all clusters except 
GC1 most likely. This constitutes the long-term timing problem in 
the sense of hypothesis |(I)lfrom page|2] which is only really a prob- 
lem if one assumes that the present globular cluster distribution of 
Fornax is representative of the distant past. 

Only our large-core model (LC) avoids this problem com- 
pletely. Rather than dynamical friction, this model shows 'dynami- 
cal buoyancy', when clusters are pushed out of the core. Following 
lTremaine & Weinberg dl984l) . this is likely caused by this model's 
inverted distribution function, when d//d£ > over a significant 



range of orbital energies E. Such m odelfl are generally thought 
to be unstable dBinnev & Tremaind 120081 §5.5). Hence, it seems 
unlikely that such a model emerges from any reasonable formation 
mechanism. However, our A'-body models appears to be stable, and 
some further research is needed to clarify the cause for 'dynami- 
cal buoyancy', i ts relation to the core stalling effect (reported by 
iRead et alj 2006). and its prospects in natural systems. 

For cusped mass models (IC & SC), such as predicted by 
CDM cosmogony via dissipationless simulations, clusters GC3 or 
GC4 will sink into the centre of Fornax within 1-2 Gyr with ~ 90% 
probability (in the sense that solutions where this does not occur 
cover only ~ 10% of the possible orbital phases and projections). In 
fact, we estimate the probability that no cluster obtains r apo < lOOpc 
within 2 Gyr to be no more than 2% in Table[3] This is largely inde- 
pendent of the assumed globular-cluster orbital distribution func- 
tion and constitues the more severe immediate timing problem in 
the sense of hypothesis |(ii)| from page [2] It implies that ;/ Fornax 
indeed has a cusped density profile, our cosmic epoch of observa- 
tion is necessarily very special. 

For a shallow cusp model (WC), only the most massive cluster 
GC3 has a 90% probability to reach, within 1-2 Gyr, r apo < 200pc 
significantly less than its current projected radius of R p = 430pc. 
For this model, no cluster can reach r apo < lOOpc within 2 Gyr and 
the chances that no cluster will reach r apo < 200pc within 1 or 2 Gyr 
are, respectively, 92% and 32% (see Table [3}- These numbers are 
not unlikely and hence avoid the long-term timing problem. 

6.2.1 A steady-state solution 

Our finding thus suggest two possible solutions to the timing prob- 
lem. First, Fornax has a large core (perhaps between models LC 
and WC) and dynamical friction is slow or has stalled a long time 
ago. In this case, Fornax may have been on its current orbit for a 
Hubble time with its globular cluster system hardly evolving. In 
this case, the consistency of the cluster distribution with the stel- 
lar distribution (as discussed in §|2.2| l cannot be a coincidence, but 
hints at a common formation scenario. 

6.2.2 An evolving solution 

In the second solution Fornax has a small core or shallow cusp (as 
model WC) and dynamical friction is still ongoing, albeit slowly 
enough that the absence of a central nucleus in Fornax (or in fact 
any other dSph) is perfectly plausible. In this case, the clusters must 
have been further away from Fornax in the past than today, and the 
current (weak) consistency of their distribution with that of the stars 
is just a coincidence. Also, a Hubble time ago the clusters most 
likely were more than the current tidal radius of 1.8-2.8 kpc away 
from Fornax. This in turn requires Fornax did not orbit the Milky 
Way for a Hubble time on its present orbit. However, even a simple 
adiabatic evolution of Fornax's orbit may be sufficient to solve this 
problem. For example, for a slowly growing Galaxy with always 
flat rotation curve the peri-centric tidal radius of Fornax evolves 
like r... oc i/~.„ (even when neglecting mass-loss of Fornax due to 

tid circ v to to 

Galactic tides). 

This second solution appears more natural and also fits with 
the weak indication of mass-segregation, as would be induced by 

4 For self-gravitating systems, an inverted distribution function typically 
occurs if the transition between a near-constant density core and steep 
power-law decay is fast, i.e. when rj in equation {5) is larger than ~ 2. 
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dynamical friction, in the current mass-radius relation (see Fig.[T]l. 
However, this model implies that the globular clusters have not 
formed within Fornax, but are most likely accreted. One may, of 
course, consider these two solutions as the extreme ends of a conti- 
nuity of solutions with various degrees of cusp strengths and hence 
dynamical friction effects. 



6.3 The case of GC1 

An interesting aspect relates to the cluster GC1. 
Penanubi aTwalker & Gilmorel d2009l) demonstrated that, uniquely 
of all Fornax's globular clusters, GC1 would be tidally disrupted 
if it fell to the centre of Fornax. While we find that GC1 does not 
sink to the centre of Fornax in almost all of our models, this still 
leaves us with a puzzle. Why should the one cluster vulnerable to 
tides be on an orbit where it would hardly ever suffer disruption? 
For our steady-state solution to the timing problem above, this 
puzzle can be resolved by the postulation that Fornax once had a 
richer globular-cluster system and we only see the survivors. Such 
survivors are either massive enough or on remote orbits to avoid 
tidal disruption. However, a high (initial) frequency of globular 
clusters (of ~ 10 4 ~ 5 M o ) appears rather implausible for a small 
galaxy like Fornax. 

For the evolving solution to the timing problem, on the other 
hand, the fact that GC1 would be disrupted poses no problem at 
all. In this picture, low-mass clusters, such as GC1, would not be 
dragged down much, and there is no need to postulate a large early 
population of clusters. 
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APPENDIX A: NUMERICAL CONVERGENCE 

In order to ensure that our simulations do not suffer from numeri- 
cal noise we ran simulations with two different mass models, one 
cusped and one cored; each with two different orbits, one circular 
and one eccentric; and each of these at three different resolutions: 
with N = 4x 10 5 , 10 6 and 4x 10 6 . We then compared the evolution 
in each case of one cluster over 10 Gyr. 

The evolution of the orbital radius of a single cluster mov- 
ing on an eccentric orbit in model WC (see table [2} is shown in 
Fig. lAll It can be seen that orbital evolution is very similar for all 
resolutions. In particular, the decay of the orbit follows the same 
timescale, with the time and radius of the first-stalling of the clus- 
ter being the same. It has been shown that the two-body noise in a 
simulation can cause the cluster orbit t o precess and c ause artificial 
decay of the orbit once in the core dRead et all 2006). We selected 
the a bove combination of orbit and density profile precisely be- 
cause iRead et all ( 120061) showed that convergence is most difficult 
for an eccentric orbit in a cored halo. This is the case where nu- 
merical friction caused by orbit precession has the largest effect on 
orbital decay. The simulations shown above give a strong indica- 
tion that such effects are not significant at even lower resolutions 
than the one used for the main body of this work. Our simulations 
are well converged. 



APPENDIX B: THE PROBABILITY FOR R^R P 

The fraction of orbital phases and projections for which R^R p on 
a given orbit is identical to the probability for R ^R p , when R is 
computed for that orbit at random orbital phase and projection. 

The probability for a cluster with orbital energy E and angular 
momentum L to be found at a radius between r and r + dr is 



p(r\E,L) dr = 



2dr 



T r v r 

and zero otherwise. r per i 
v 2 = 2[E-1>(r)]-L 2 /r 2 



if ^ peri ^ f ^= r apo 



(Bl) 



. are the roots of the radial velocity 

(B2) 



and the radial period is given by 
dr 



7V = 2 



p P o dr 



(B3) 



The probability that a cluster at radius r has projected radius 
between R and R + dRis 



p(R\r)dR- 



RdR 



if R<r 



(B4) 




Figure Al. The evolution of the radius for a single cluster moving on an 
eccentric orbit in a halo with a near-cored density profile (model WC from 
Table|5J, but realised with different particle numbers as indicated. Some of 
the deviation at later times is due to the uncertainties in determining the 
centre of the A'-body model. 

and zero otherwis^f]. Thus, the probability for a cluster with given 
orbit to have projected radius between R and R + dR is 



p(R\E,L)dR = dR 



p(R\r)p(r\E,L) dr 



(B5) 



for R < r apo and zero otherwise. From this, we can work out the 
finite probability that a cluster with given orbit has projected radius 
not greater than observed: 



p(R^R p \E,L) 



p(R\E,L)dR 

1 r JO Jim 

2 pi» dr r' 
TrJr^ VrJo 



mIK,^} r^r 2 -R 2 V r 

rain(r.R„) R ^ 



rir 2 -R 2 



(B6) 



with (•)+ = maxjO, •). One can easily verify that p(R < 0\E,L) = 
and p(R < r apo |£, L) = 1, as required. 



5 This follows from the probability density p(8)d8 = sinfffor the polar angle 
Be [0,n] and p{R\r)AR = p(8){d8l 6R)AR with R = rsind. 
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